library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # TS padding
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz

library(forecast) # TS library
library(TTR) # SMA function
library(tseries) # adf.test

# supaya semua plot memiliki theme_minimal()
theme_set(theme_minimal())

1 Dive Deeper

Gunakan data souvenir untuk melakukan fitting model dan forecasting menggunakan metode exponential smoothing. Data tersebut merupakan data penjualan souvenir per bulan selama 7 tahun dari bulan Januari 1987 sampai Desember 1993.

1.1 Import Data

Load data menggunakan fungsi scan() untuk emmbaca data menjadi sebuah vector, lalu lakukan visualisasi dengan base plot()

# import data
souvenir <- scan("data_input/fancy.dat")

# line plot
plot(souvenir, type = "l")

1.2 Time Series Object

Buatlah object Time Series ts(data, start, frequency), lalu visualisasikan dengan autoplot()

# object TS
souvenir_ts <- ts(data = souvenir,
                  start = c(1987, 1),
                  frequency = 12)

# autoplot
souvenir_ts %>% 
  autoplot()

Insight: Apakah TS termasuk tipe additive/multiplicative? multiplicative

1.3 Exploratory Data

Lakukan decomposition sesuai dengan tipe additive/multiplicativenya, lalu buatlah bar plot untuk menganalisa komponen seasonalitynya.

# decomposition
souvenir_dc <- souvenir_ts %>%
  decompose(type = "multiplicative") 

souvenir_dc %>% 
  autoplot()

# seasonality analysis
data.frame(
  Month = factor(month.abb, levels=month.abb),
  Seasonality = souvenir_dc$seasonal
) %>% 
  distinct() %>% 
  ggplot(aes(x = Month, y = Seasonality)) +
  geom_col()

Insight: Bagaimana trend dan seasonality dari penjualan souvenir?

  • Trend: cenderung naik
  • Seasonality: penjualan memuncak di akhir tahun, yaitu Nov dan Dec

1.4 Cross Validation

Lakukan train-test splitting, dimana 1 tahun untuk data test dan sisanya untuk data train

# data test
souvenir_test <- tail(souvenir_ts, 12)

# data train
souvenir_train <- head(souvenir_ts, -12)

1.5 Fitting Model

Buatlah minimal 2 model exponential smoothing dengan fungsi berikut:

  • HoltWinters()
  • ets()

Note, pemilihan model berdasarkan:

  • ada atau tidaknya trend, seasonality
  • tipe additive/multiplicative
souvenir_hw <- HoltWinters(souvenir_train, seasonal = "multiplicative")
souvenir_hw_add <- HoltWinters(souvenir_train %>% log(), seasonal = "additive")
souvenir_ets <- ets(souvenir_train, model = "ZZZ")

1.6 Forecasting

Gunakan forecast() untuk masing-masing model yang Anda buat sebelumnya:

souvenir_hw_f <- forecast(souvenir_hw, h = 12)
souvenir_hw_add_f <- forecast(souvenir_hw_add, h = 12)
souvenir_ets_f <- forecast(souvenir_ets, h = 12)

Visualisasikan data actual (train dan test) dengan hasil forecasting

p <- souvenir_train %>% 
  autoplot(series = "Train") +
  autolayer(souvenir_test, series = "Test") +
  autolayer(souvenir_hw_f$mean, series = "HW Multiplicative") +
  autolayer(souvenir_hw_add_f$mean %>% exp(), series = "HW Additive") +
  autolayer(souvenir_ets_f$mean, series = "ETS")

plotly::ggplotly(p)

1.7 Model Evaluation

accuracy(souvenir_hw_f, souvenir_test)
#>                      ME     RMSE      MAE         MPE     MAPE      MASE
#> Training set   189.1763 2200.411 1512.406  -0.1299785 13.13863 0.4038682
#> Test set     -4060.1992 4541.688 4060.199 -21.1511565 21.15116 1.0842231
#>                   ACF1 Theil's U
#> Training set 0.1829319        NA
#> Test set     0.5148249 0.6278038
accuracy(souvenir_ets_f, souvenir_test)
#>                       ME     RMSE      MAE         MPE     MAPE      MASE
#> Training set    84.47865 1805.020 1228.946  -0.5220307 11.80610 0.3281741
#> Test set     -1019.41169 3790.485 3389.888 -10.2903652 18.22484 0.9052254
#>                   ACF1 Theil's U
#> Training set 0.2111953        NA
#> Test set     0.7158307 0.5309041
# cara mengevaluasi apabila data multiplicative diperlakukan sebagai additive
accuracy(souvenir_hw_add_f$fitted %>% exp(), souvenir_train) # training
#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -55.25421 2302.798 1552.995 -1.715425 13.25152 0.1180214 0.3660692
accuracy(souvenir_hw_add_f$mean %>% exp(), souvenir_test) # testing
#>                 ME     RMSE      MAE      MPE    MAPE        ACF1 Theil's U
#> Test set -5769.926 6138.801 5769.926 -25.9165 25.9165 -0.04706218 0.6797808

Kesimpulan: Berdasarkan nilai error di atas, model yang terbaik berdasarkan metric MAPE adalah souvenir_ets yaitu ETS(M,A,M).

1.8 Next Step

Bagian ini opsional, namun penting untuk diketahui

Misalkan selanjutnya ada kebutuhan untuk melakukan forecasting periode Jan 1994 - Dec 1994 (satu tahun kedepan dari data souvenir). Maka practice yang baik, biasa kita re-train semua data kita menggunakan best model pada tahap sebelumnya

# misal pada tahap sebelumnya, model terbaik adalah dengan fungsi ets() model ZZZ
souvenir_ets_full <- ets(souvenir_ts, model = "ZZZ")
souvenir_ets_full_f <- forecast(souvenir_ets_full, h = 12)

# visualisasi
souvenir_ts %>% 
  autoplot(series = "Full Data") +
  autolayer(souvenir_ets_full_f, series = "Forecast 1 year")

Model souvenir_ets hanya di train dari Januari 1987 sampai Desember 1992, sedangkan souvenir_ets_full di train dengan keseluruhan data (sampai Desember 1993).